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A NUMERICAL METHOD OF DETECTING SINGULARITY OF A MATRIX 


M. La Porte and J. Vignes 


I. Introduction /7 3* 

From a theoretical point of view, a matrix is singular if, 
and only if, the determinant of its coefficient is zero. Clas- 
sical numerical algorithms permit calculating the value of a 
determinant. One such calculation carried out on a computer does 
not generally give a zero value for the determinant when the 
matrix is singular. In effect, the result obtained is vitiated 
due to the limited arithmetic precision of a computer. Then, 
there is always disagreement between the analytical result and 
the numerical result obtained on the computer. From this fact, 
it is apparent : 

— a matrix, analytically non-singular, can appear as singular, 

— a matrix, analytically singular, can appear as non-singular. 

It is easy to give examples illustrating this disagreement. 

Let us consider the matrices: 

^ + i|' ;|- 

The matrix A is regular and its determinant has a value of 
of 10"'^. The B matrix is singular. 

Let us treat these two matrices on an imaginary computer 
operating to 7 significant points, arithmetically at a normalized 
floating comma with truncation. 

The matrix A is presented in the form: 


^♦Numbers in the margin indicate pagination in the foreign text. 
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all of the being equal to 1. Then one has A ^ A. The value 
of determinant A, calculated with any algorithm, is zex’o. This 
matrix appears then as singular when it is not analytically. 

The B matrix is presented in the form B = 5 since its / 7 ^ 

coefficients fall correctly in the machine. Its determinant, 
calculated using a hausslan algorithm with study of the maximal 
center-point by column, has a value of: 

/I* =3 X io-«; 

Then, this matrix appears as non-singular when it is 
obviously singular. 

These two examples show the difficulty in detection, by an 
exclusively numerical method, of the singularity of one matrix. 

These studies bearing on the propagation of numerical errors 
during calculation and their incidence in linear algebra have 
been made particularly by: Forsythe [1], Gastinel [2], Golub [3], 

Householder [4, 5], Wilkinson [7-10]. 

II. Definition of Numerical Singularity 

If the calculation of a determinant of an analytical 
singular matrix has a value which is not zero, one can confirm 
that this value is only the result of a cumulative effect of 
errors caused by the arithmetic of the computer. From this fact, 
this value is not significant and, consequently, it must be con- 
sidered as mathematically zero. 

Moreover, if the calculation of the determinant of an 
analytically non-singular matrix gives a zero value, one can say 
that on this computer this matrix appears as singular whereas it 
does not analytically. This necessitates going beyond the 
concept of analytical singularity in order to arrive at the 
concept of numerical singularity; for this we propose the following 
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definition; 


Definition. A matrix is numerically singular when the value 
of its determinant, calculated on the computer, is: 

-‘-either zero, 

--or not significant. 


Let us consider a matrix A in which the coefficients 
the algebraic values: 



are 


*«/ = [«</]. 


( 1 ) 


the 


“ij 

6 ]. 


Let us designate as A the image on the machine A and by 
images of . In general, the A^j are not strictly equal to 
due to error engendered by the periodic block operator [1, 5, 


They may be: 

— Det(A) the algebraic value of the determinant of A, 

— Aj the value calculated on the computer of determinant A 
— e the absolute error committed on Aj 

If the error e is of the same order of magnitude as the 
value Aj, one can confirm that this latter is not significant. 
Consequently, the matrix must be considered as numerically singula'^ 

If e is smaller than Aj, then the matrix is not singular. 

It is impossible to calculate e by (2) since the value of 
Det(A) is unattainable. Here we will present a method which, by 
statistical. methods, makes it possible to estimate e. 

III. Origins of Error 


/75 
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Absolute error e, given by equality (2), has two distinct 
origins. 

a) error in calculation. This corresponds to the cumulative 

effect of errors caused by each of the elementary operations 
carried out on the computer during calculation of [1, 6, 7]. 

b) error in the coefficients. The coefficients of matrix 
A are not in general equal to the coefficient a^j of the algebraic 
matrix A since they depend on part of the operator periodic block, 
and the other part on the orig’*^ of • Three cases can be pre- 
sented : 

Case 1 . The a^j are known exactly and go correctly into the 
machine. In this case; A = A. 


Case 2 . The a^^ are known exactly but do not fall correctly 
in the machine. A is an approximation of A, the error on the 
coefficient A^j proceed only from a single operator periodic 
block . 


Case 3 . The a^j are known with a certain error; this is the 
case, for example, when the a^j result from experimental measurment . 


In this case, the periodic block error is negligible and the A 


are affected by the same error as the a 


IJ 


1J‘ 


We are going to determine, by a statistical method, a mean 
estimate of error e, starting with the populations created 
differently following the origin of error. The elements of these 
populations are the different numerical values of the the same 
determinant . 


— if the a^j correspond in the first 
calculation has to be considered and the 
is obtained by a "permutation" method. 


case, only error in 
corresponding population 


— If the correspond to the second and third cases, the 

error in the coefficients and the error in the computation must 
be taken into account and the corresponding populations obtained 
by a "permutation-perturbation" method. 


IV. Permutation Method (Theoretical Aspect ) 

This method applies to the first case, that is to say, when 
the coefficients of A are the exact image of coefficients 
a^^j and A. 


( 3 ) 

During calculation on a computer of the value of Aj, the absolute 
error committed e is uniquely caused by the limited arithmetic 
precision of the machine. Due to this, it depends on the order 
in which the operations are carried out. 

Also, in permuting the columns of A, one changes the order 
of operations and if one calculates the determinant, one obtains 
different value of Aj. 

By making all of the possible permutations of the column of / 7 6 
matrix A in n order, one obtains n! matrices. Machine compu- 
tation of n\ corresponding determinants produces a population 

Dj. 

Z)| = • • • » • • •)• 

All of the values of A^ are also represented by one or 
another value of Det(i4). 

Also, the permutation method makes it possible to theoretically 
obtain a population Dj of the cardinal: 

CardD|=»l ... 
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V. Permutation-perturbation Method (Theoretical Aspect ) 

This method is applied to the second and third cases, that 
is to say, when the coefficients a . do not fall correctly in 
the machine. Then one has: 




( 5 ) 


When one is in case No. 2, the correspond to periodic 
block error and belong to the P(a) population defined in [6] 

§ 2.2. Then, one has: 


+«</) 


( 6 ) 


When one is in case No. 3» the X^j correspond to those 
experimental errors in which one assumes a higher boundary e 
as knov.n. Then, one has: 


ij 




(7) 


In the case defined by (6), only two representations on the 
computer exist for each of the values of a.,, the one by 

_ 4 - ^ J 

efficiency, , the other by excess and one has: 




( 8 ) 


If one works out all of the resulting matrices from a com- 
bination of the two possible states of each of the coefficients, 
one obtains 2^ matrices. By applying the permutation method to 
each of these matrices, one obtains a population D 2 : 




Card • 

In the case defined by (7), for each of the a^j , there 
exist Pj^j representations on the computer. 


( 9 ) 
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If one works out all of the matrices resulting from a 
combination of possible states from each of the coefficients, 
one obtains N matrices: 


N-// riPn^ 

<-i 


( 10 ) 


In applying the computation method to each of these matrices, 
one obtains a population D3: 

. 


Card 1 17 ■ 17 ^ 1 #* . 

I tMl 


( 11 ) 


VI, Evaluation of t (Theoretical Aspect) 

Using the elements of one of the populations we find 
above : 


^h) ^^ 2 ) 

we will get a mean estimate of e of absolute error e defined by: 

with Aj = the caclulated value on the computer of determinant A 
Det (A) = the exact value of the determinant of A. 

Here we can make the following hypothesis: 

Hypothesis . During calculation of A^ on the computer, the 
calculation errors of data have behaved indifferently in one 

I' 

direction or the other. Consequently, the value of Det(A) can 
be considered as an element of any of the D population. 

We will call t and 6^, respectively, the mean and the 
variants of this population. Under the preceding hypothesis, 
the mean of error square e is expressed by: 




( 1 ^) 


By definition, we suppose: 




( 15 ) 


or 




( 16 ) 


The value of t compared to that of Aj gives us information 
on the conditioning of the matrix. In the measurement where Aj 
and e are not zero, one can always write: 


I 

TjT 


= 10 -® 


( 17 ) 


C represents the exact number of significant figures of A^. Of 
course, the value of C can never exceed the number C „ of sig- 
nif leant decimals in the arithmetic of the computer. If p 
designates the number of bits of the mantissa representing a 
normalized floating binary comma, C is defined by: 

niaX 


2 * = 10 ®—. 


( 18 ) 


Then, three cases can be considered: 

— C is in the neighborhood of C : the matrix is 

max 

numerically well conditioned. 

— C is in the neighborhood of zero: the matrix is 

numerically singular. 

— Beteen the two ext. “’erne cases, the matrix without being 
singular is not numerically well conditioned. 
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We will say that the matrix Is numerically singular If 
the value of Its determinant Is known with less than an exact 
significant figure, that Is to say. If C < 1. 

VII. Practical Aspect of the Permutation and Perturbation Methods 

In practice. It Is Impossible to work out the totality of 
the population elements previously defined. We propose an 
algorithm which, with a number of limitations, makes It possible 
nevertheless to determine the number C. 

Using any algorithm from calculation of the determinant, 
one defines the value of Aj from the determinant of matrix A and 
the value of A 2 of the determinant of the matrix deduced from A 
by central symmetry. The Interesting thing about this choice 
Is the fact that during calculation of Aj and A 2 , the errors are 
propagated In a very different way. 

Using these two elements, one calculates by (17) the value 

of C: 


— If C < 1, then the matrix Is numerically singular, 

— If C > 1, the algorithm must be pursued by creating new 
elements of the D population until one obtains a stationary con- 
dition for the entire part of C. The successive elements of D 
are obtained by random permutation of the columns of matrix A 

and, should the occasion arise, by random perturbation of the 
A^j coefficients. This perturbation consists of: 

— In case No. 2, by random adjustment of a zero or a 1 at 
the last bit of the mantissa, 

— In case No. 3, by replacing A^^j by: 




( 19 ) 


randonily takes the value -1 or +1. 

VIII. Agreeement Between Theory and Practice 

The permutation-perturbation method has been applied to a 
very large number of matrices of all orders and all conditions. 
This method has never been deficient. We will content ourselves 
with giving three important examples here. 

VIII-1. Application of the Hilbert Matrix 

The Hilbert matrix is a symmetrical matrix in which the terms 
are the inverse of successive integrals. 




i -1/2 . i/5 

i/2 i/3 i/4 

i/n i/(n + i) 


Hn 

i/(» + i) 
i/(2*-i) 


Using the Gaussian algorithm, one calculated on a CDC 7600 
computer for n variants of 2 to 13, the value of then the value 
of C by the permutation and perturbation method (case No. 2). 

Likewise, one knows the algebraic expression of the deter- 
minant of A : 
n 


with 


Det(jr.) 




V(")“ // (»i). 




( 20 ) 

one can deduce from this the exact value of C, or C*, defined by: 

(21) 


D®t(^r>) “ ii| 4n— c* 
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The results presented in Table i show that: 


— on the one hand, perfect agreement between theory and 
practice (with values of C and C* are similar), 


— on the other hand, the small number of calculations of 
determinants necessary for obtaining the results. 


Table I 


• 

Del(jrj 


Number of 
calculations 

Either nart ** 
c of c* 

1 

8.J X lO"* . 

, 8.3 X IO-* 

■ 4 

1 » . 

13 

14 . 1 1 

3 

4.6 X ior* , 

4.6 X IO-* 

. 3 

1 . J 

12 

12 ■ 

4 

l.6xl(r’ 

1.6 X 10"’ 

3 


11 

, <0 . 

5 

3.7 X l0-‘» 

3-7 X I0'“ 

3 

* * * 

9 


6 

$ 3 X l(T« 

5 3 X 10-“ 

3 

■ 

8 

. • g • . f • 

7 

4.8 X 10“* 

4.8 X I0-" 

4 


7 

7 

1 

2 7 X IO-“ 

2.7 X IO-“ 

4 


5 

6 

9 

9 7 X to-** 

9 7 X IO-“ 

3 


4 

4 

10 

2.1 X l(T“ 

2.1 X IO-“ 

3 


3 

2 

II 

3.0 X IO-* 

2.3 X IO-“ 

3 


1 

• 

12 

2.6 X 10-" 

••4.5 X IO-“ 

2 


0 

0 5 

tj 

1.4 X 10-** 

-7.8 X IO-" 

3 


0 

0 5 


S: numerically singular matrix. 

OiUijlN’AI, 

VIII-2. Application of Least Squares to a Matrix (jUALITlfi 


Let us consider the matrix 








V. 


s. 


with 


which one encounters in certain polynomial adjustments by 
least squares. 


The results obtained for N = 20 and p variations of 1 to 12 / 80 

are presented in Table II. This shows the good agreement between 
C and C* and also shows that a matrix can be numerically singular 
Itself if the order of magnitude of its determinant is very high. 
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Table II 




Number of A 


5* 



• * 


calculations 

c 

c* ■ ‘ : 

1 

1.6x |y| 

1.6 X 10* 


13 

1 12 : 1 r. r- 


a • 


3.6 X 10* 

'■ '*• 4 ?l| Il'.M. .|l 1 1 , 

13 *. 

ia ’■ t- I - 1 

f| . ! . , 

1 

4 

IlKloX 
31 X 10» 

3.8 X 10*1 

' 4 * *» . 1 

i 

!l 

$ • 

1.7X10** ‘ 

‘ 1.7x10** 

' ' 3 It I, 

10 

1. 1. • 


6 

1.9 X10** • 

.• 1.9 X 10** 

- *• 3 .. ,.i 


8 1 » ' 

7 

S.a X io»* 

sa X 10** 


7 • 

7 . •'*» 

8 

9 

3.4 X IC^« • 
5.1 X 10’* 

3.4 X 104* 
5.1 X 10** 

^ ‘ • •• 1 , 

3 

5 t 

4 

5 

10 

1.6X10** 

1.5x10** 

3- 

3 

a 

11 

1.0 X 10*** 

-3.6 X 10*** ' 

a . ■ 1 

1 

A 1 

1 

11 

i.axio*** 

-1.8 X 10*** 

a 

0 

0 

0 5 

0 c 


S: numerically singular matrix. 


original Pirp ro 
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VIII-3. Application to Algebraically Singular Matrices 


We have processed more than ten thousand matrices on the 
CDC 7600 computer with order of magnitude varying from 2 to 100, 
the coefficients > drawn at random being the orders of mag- 
nitude between 10“® and lO’’"^ . One has found these singular matrices 
by replacing the last line with the sum of the other lines. 


The permutation-perturbation method has always made it 
possible to draw a conclusion as to the singularity of these 
matrices. Table III shows that this result has been obtained 
using a number of elements, never exceeding 3. On the other 
hand, the higher the order, the more the mean number of necessary 
elements tends toward 2. 


IX. Conclusion 

The permutation and perturbation method has always given 
satisfaction and has never been defective and therefore one can 
designate it as certain. On the other hand, from this study 
it stands out that : 
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Table III 



— when the matrix Is not numerically singular, a population 
of 3 or 4 elements is sufficient to assure and know the numerical 
conditioning of the matrix. 

— when the matrix is numerically singular, its singularity 
is generally found from a population of two elements. 

One can then say that this method is certain and efficient. 

It can be used for numerous problems of linear algebra. Also, 
it makes it possible to control the intrinsic value X of matrix 
M verifying that the M-XI matrix is singular. 

Applied here to a problem of linear algebra, this method 
can be generalized for cases of nonlinear algebra because its 
concept basically is to execute the same algorithm many times 
propagating the errors in various ways finally to obtain results 
which vitiate the various errors. 

ORIGINAL PAGE IS 

This generalization consists of: OF POOR QUALITY” 

— on the one hand, changing the order of conducting / 8l 

elementary operations (permutation), 

— on the other hand, taking its value by deficiency or excess 
(perturbation) as the result of each elementary operation. 
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It permits introducing a new concept of computers making 
them capable of combining an estimate of their precision in the 
total results. 


Ii4 
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